// math_utils.h

#ifndef DODO_MATH_UTILS_H
#define DODO_MATH_UTILS_H
#include <cmath>
#include <algorithm>
#include <stdexcept>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <boost/tuple/tuple.hpp>

#include <quigon/basic/constants.h>
#include <quigon/basic/basic_utils.h>
namespace quigon
{

class Float_Equal
{
public:
	Float_Equal(double e0 = 1e-5) : epsilon(e0) {}

	bool operator()(double a, double b) const
	{
		if ((a - b) <= epsilon && (b - a) <= epsilon)
			return true;
		else
			return false;
	}

private:
	double epsilon;
};

const Float_Equal float_equal;
const Float_Equal float_equal1(1e-1);
const Float_Equal float_equal2(1e-2);
const Float_Equal float_equal3(1e-3);
const Float_Equal float_equal4(1e-4);
const Float_Equal float_equal5(1e-5);
const Float_Equal float_equal6(1e-6);
const Float_Equal float_equal7(1e-7);
const Float_Equal float_equal8(1e-8);
const Float_Equal float_equal9(1e-9);
const Float_Equal loosely_equal(1e-3);

// require Tx, Ty:
// operator <
// mathematic + - * /
template<class Tx, class Ty>
Ty interp(Tx x1, Ty y1, Tx x2, Ty y2, Tx x)
{
	try
	{
		using namespace std;

		if (x2 < x1)
		{
			swap(x1, x2);
			swap(y1, y2);
		}
		if (x < x1)
			return y1;
		if (x2 < x)
			return y2;

		if (x1 == x2)
		{
			if (y1 == y2)
				return y1;
			else
				throw invalid_argument("Interpolating with x1 == x2");
		}

		return (y2 - y1) * (x - x1) / (x2 - x1) + y1;
	}
	catch (...)
	{
		func_err("interp");
		throw;
	}
}

template<class Ty>
Ty interp(
        boost::posix_time::ptime x1, Ty y1,
        boost::posix_time::ptime x2, Ty y2,
        boost::posix_time::ptime x
)
{
	try
	{
		using namespace std;

		if (x2 < x1)
		{
			swap(x1, x2);
			swap(y1, y2);
		}
		if (x < x1)
			return y1;
		if (x2 < x)
			return y2;

		if (x1 == x2)
		{
			if (y1 == y2)
				return y1;
			else
				throw invalid_argument("Interpolating with x1 == x2");
		}

		return
		        (y2 - y1) *
		        (
		                (x - x1).total_nanoseconds() /
		                static_cast<double>((x2 - x1).total_nanoseconds())
		        ) + y1;
	}
	catch (...)
	{
		func_err("interp");
		throw;
	}
}

// This is non of the maths, but to keep all the interp together, I put it here
template<class Ty> Ty interp(
        const std::string& x1, Ty y1,
        const std::string& x2, Ty y2,
        const std::string& x
)
{
	throw std::invalid_argument("Interpolating with string type.");
	return y1; // hoax for the compiler...
}
template<class Tx> std::string interp(
        const Tx& x1, const std::string& y1,
        const Tx& x2, const std::string& y2,
        const Tx& x
)
{
	throw std::invalid_argument("Interpolating with string type.");
	return y1; // compiler stuff
}

// requirement: *x and *y can do arithmetics with double and results double.
template<class Iter_x, class Iter_y>
boost::tuple<double, double, double>
linear_fit(Iter_x beg_x, Iter_x end_x, Iter_y beg_y, Iter_y end_y)
{

	double x_bar = 0, y_bar = 0;
	double lxx = 0, lyy = 0, lxy = 0;
	Iter_x iter_x = beg_x;
	Iter_y iter_y = beg_y;
	size_t count = 0;
	for (; iter_x != end_x && iter_y != end_y; ++iter_x, ++iter_y)
	{
		x_bar += *iter_x;
		y_bar += *iter_y;
		++count;
	}
	x_bar /= static_cast<double>(count);
	y_bar /= static_cast<double>(count);

	for (iter_x = beg_x, iter_y = beg_y; iter_x != end_x && iter_y != end_y; ++iter_x, ++iter_y)
	{
		lxx += (*iter_x - x_bar) * (*iter_x - x_bar);
		lyy += (*iter_y - y_bar) * (*iter_y - y_bar);
		lxy += (*iter_x - x_bar) * (*iter_y - y_bar);
	}
	double rxy = lxy / sqrt(lxx * lyy);
	double b = lxy / lxx;
	double a = y_bar - b * x_bar;

	return boost::make_tuple(a, b, rxy);
}
// requirement: container has value_type, begin(), end(), size()
template<class Container_x, class Container_y>
boost::tuple<double, double, double>  //a, b, rxy
linear_fit(const Container_x &x, const Container_y &y)
{
	return linear_fit(x.begin(), x.end(), y.begin(), y.end());
}

}

#endif

